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ABSTRACT 

Hierarchical assembly models predict a population of supermassive black hole (SMBH) 
binaries. These are not resolvable by direct imaging but may be detectable via periodic 
variability (or nanohertz frequency gravitational w aves). Foll owin gpo ur detection of a 
5.2 year periodic signal in the quasar PG 1302-102 (Graham et al. 2015), we present a 
novel analysis of the optical variability of 243,500 known spectroscopically confirmed 
quasars using data from the Catalina Real-time Transient Survey (CRTS) to look for 
close (< 0.1 pc) SMBH systems. Looking for a strong Keplerian periodic signal with 
at least 1.5 cycles over a baseline of nine years, we find a sample of 111 candidate 
objects. This is in conservative agreement with theoretical predictions from models of 
binary SMBH populations. Simulated data sets, assuming stochastic variability, also 
produce no equivalent candidates implying a low likelihood of spurious detections. The 
periodicity seen is likely attributable to either jet precession, warped accretion disks 
or periodic accretion associated with a close SMBH binary system. We also consider 
how other SMBH binary candidates in the literature appear in CRTS data and show 
that none of these are equivalent to the identified objects. Finally, the distribution 
of objects found is consistent with that expected from a gravitational wave-driven 
population. This implies that circumbinary gas is present at small orbital radii and is 
being perturbed by the black holes. None of the sources is expected to merge within at 
least the next century. This study opens a new unique window to study a population 
of close SMBH binaries that must exist according to our current understanding of 
galaxy and SMBH evolution. 

Key words: methods: data analysis — quasars: general — quasars: supermassive 
black holes — techniques: photometric — surveys 


1 INTRODUCTION 


Supermassive black hole (SMBH) binary systems are an 
expected consequence of hierarchical models of galaxy for¬ 


mation (Haehnelt & Kaufmann 2002 Volonteri, Madau & 


Haardt 2003) and an important sources of nanohertz fre¬ 


quency gravitational waves. Theoretically, they are seen to 
evolve through three stages (for a recent review, see Colpi 
(2014) and references therein): in the first, the SMBH pair 
sinks towards the center of the newly merged system via dy¬ 
namical friction (the SMBHs feel the collective effect of the 
stellar distribution). The binary continues to decay due to 
interactions with stars in the nuclear region on intersecting 


* E-mail:mjg@caltech.edu 


orbits with the binary. These stars carry away energy and 
angular momentum. The merger process may stall in this 
phase (at separations of 0.01 - 1 pc) owing to the depletion 
of stars in the nuclear region. However, other factors such 
as the binary mass, the presence of cold gas in the nuclear 
region, or deviations from axial symmetry in the galaxy rem¬ 
nant can affect the length of the potential stall. Depending 
on these factors as well as the stellar distribution, the binary 
will generally spend most of its lifetime in this stage. Finally, 
if the binary separation decreases sufficiently through these 
processes, the remaining orbital angular momentum of the 
pair is efficiently dissipated by the emission of gravitational 
radiation, leading to an inevitable merger. At coalescence, 
the merged SMBH receives a kick velocity and may recoil 
and oscillate about the core of the galaxy or even be ejected. 
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SMBH binary pairs have been detected at kilopar- 
sec separations, e.g., NGC 6240 at a separation of 1.4 
kpc by Chandra ( Komossa et al.||2003| , but the obser¬ 
vational evidence for close (subparsec) pairs is tentative. 
The highest milliarcsecond angular resolution imaging with 
VLBA/VLBI radio observations is only feasible for the very 
nearest systems - a 1 pc separation pair at a distance of 100 
Mpc (z ~ 0.024) subtends 2 milliarcsecs. Rodriguez et al. 


sion lines from the two BLRs associated with the SBMHSs 
will be resolved since their separation is likely to be smaller 
than their widths. The optical and UV spectral features are 
therefore no longer directly related to the period of the bi¬ 


nary (Shen V Loeb 2010). However, variability related to 
the periodicity of the accretion flows onto the binary may be 
directly detectable (Artymowicz V Lubow 1996; Hayasaki, 


(2006) report the discovery of two resolved compact, vari¬ 


able, flat-spectrum radio sources with a projected separation 
of 7.3 pc in the radio galaxy 0402+379 (z = 0.055). The best 
known close SMBH binary candidate, OJ 287 (Valtonen et 
al. 2008), shows a pair of outburst peaks in its optical light 
curve every 12.2 years for at least the last century, which 
has been interpreted as evidence for a secondary SMBH per¬ 
turbing the accretion disk of the primary SBMH at regular 
intervals ( |Valtonen et al.||2011 ). Most other claims of pho- 


|Mineshige &; Ho||2008). 

As noted earlier, there have been reports of 
(quasi-)periodic behavior detected in the long-term multi¬ 
wavelength monitoring of particular blazars, e.g., PKS 2155- 


304, AO 0235+164 (Rani, Wiita V Gupta 


2009 


Wang 2014). 


tometric (quasi-)periodicity, e.g., Fan et al. (2002) 


Rieger 


&; Mannheim| (|2000| ; |De Paolis, Ingrosso &; Nucita (2002); 


L iu et al.J p014|), however, tend to have far shorter temporal 


baselines and do not show the same level of regularity. 

To date the search for SMBH binary systems on 
(sub-)parsec scales has therefore focused on detection by 
spectroscopic discovery, almost solely within the SDSS data 
set. It should be noted, though, that there is no unique spec¬ 
tral characteristic signature of a SMBH binary, despite sev- 


These can be separated into short timescale high-energy (X- 
ray, gamma ray) phenomena, typically with periods of a few 
hours or days, and longer timescale radio and optical be¬ 
haviors with timescales of a fraction of a year to decades. 
Given the nature of these objects, though, these are gener¬ 
ally explained as related to jet-disk interactions around a 
single SMBH rather than a SMBH binary. 

The expected period for a pair of 10 8 M© black holes 
at small separations (~ 0.01 pc) is ~10 years and there are 
now a number of synoptic sky surveys with sufficient tem¬ 
poral baselines and sky coverage for a large scale search for 
such objects based on their photometry (Pojmanski 2002 
Udalski et al.|[T993| Rau et al.||2009| Sesar et al.||2011 ). In 


fact, Haiman, Kocsis & Menou (2009) proposed systematic 


eral theoretical efforts to identify one, e.g., Shen V Loeb 

(2010 

); Montuori et al. (2011). Initial searches ( 

Bogdanovic 

et al. 

20081; Dotti et al. 2009 

Boroson & Lauer|2009 

) looked 


for double broad-line emission systems - thought to orig¬ 
inate in gas associated with the two SMBHs, where the 
velocity separation between the two emission line systems 


searches for periodic quasars in large optical or X-ray sur¬ 
veys and calculated the likely detection rates for various ob¬ 
serving strategies. The expectation is that such objects are 
rare - Volonteri, Miller & Dotti (2009) estimate that there 


should be ~10 with z < 0.7. 

In this work, we present a systematic search for periodic 


line systems with a systematic velocity offset of the broad sient Survey (CRTq 1 

Drake et al. 2009 Mahabal et al. 2011 

component from the narrow component (e.g., Tsalmantza 

Djorgovski et al. 2012 

). This is the largest open (publicly ac- 


2011 Tsai et al. 2013), under the assumption that only one 


SMBH is active and Keplerian rotation about this produces 
the shift. However, it is possible to explain such phenomena 
via a disk emitter model with a single SMBH rather than 
requiring a binary system. 

More recent searches have been based on multi-epoch 
spectroscopy, looking for temporal changes in the velocity 
shifts consistent with binary motion (acceleration effects). 
This can either be from dedicated follow-up spectroscopy 

Eracleous et al. 2012| [Tt 


of initial epoch candidates (e.g., 


carli et al. 2013) or more general searches using emerging 
collections of multi-epoch data, e.g., Shen et al, (2013); Liu 


et al. (2014); Ju et al. (2013). With typical velocity offsets 
of ~ 1000 km s -1 , searches are typically sensitive to sys¬ 
tems with mass M ~ 10 8 M© separated by ~ 0.1 pc and 
with orbital periods of ~ 300 years. They will therefore only 
capture a fraction of an orbital cycle and other astrophys- 
ical processes, such as an orbiting hotspot produced by a 
local instability in the broad line region (BLR) of a galaxy 
or outflows associated with the accretion disk, may still be 
responsible for any effect seen (Bogdanovic 2015; Barth et 
aLp0l5 ). 

At small separations (~ 0.01 pc), the size of any BLR 
is significantly larger than the semimajor axis of the binary. 
The velocity dispersion of the BLR gas bound to a single 
SMBH is higher than the velocity difference between the two 
SMBHs in a binary. It is very unlikely, therefore, that emis- 


~ 33000 deg 2 between —75° < Dec < 70° (but avoiding re¬ 
gions within ~ 10° — 15° of the Galactic plane) to a depth 
of V ~ 19 to 21.5. Time series exisl[^]for approximately 500 
million objects with an average of ~250 observations over a 
9-year baseline. 

This paper is structured as follows: in section 2, we 
present the selection techniques for identifying periodic 
quasars and in section 3, the data sets we have applied them 
to. We discuss our results in section 4. We consider the CRTS 
time series of other periodic candidates in the literature in 
section 5. We discuss our findings in section 6 and present 
our conclusions in section 7. We assume a standard WMAP 
9-year cosmology (fbv — 0 .728, 9 m — 0.272, H 0 = 70.4 


km s Mpc ; Jarosik et al. 2011) and our magnitudes are 


approximately on the Vega system. 


2 IDENTIFYING PERIODIC QUASARS 

The CRTS time series are typically noisy, gappy and irregu¬ 
larly sampled and this precludes any standard Fourier-based 
analysis of periodicity which assumes regular sampling. Vari¬ 
ant techniques that seek to regularize the data, such as 


1 http: / / crts.caltech.edu 

2 http://www.catalinadata.data 
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Lomb-Scargle, may be ~ 30% ineffective for this type of 
data ( Graham et al.||2013 ) and so are not suitable for this 
analysis. Instead we have developed a joint wavelet and au¬ 
tocorrelation function based approach to select our periodic 
quasar candidates. 


2.1 Wavelet analysis 

Wavelets are an increasingly popular tool in analyses of 
time series and are particularly attractive since they allow 
both localized time and frequency analysis. In particular, the 
power spectrum of a time series can be evaluated as a func¬ 
tion of time and the time evolution of parameters associated 
with possible (quasi-)periodic behavior determined, i.e., pe¬ 
riod, amplitude and phase. Although conventional wavelet 
analysis via the discrete wavelet transform requires regularly 
sampled data, a number of techniques have been developed 
to deal with irregularly sampled data. 

The discrete wavelet transform (DWT) of a time se¬ 
ries, At, is a convolution with the complex conjugate of the 
particular wavelet function, /, being used: 

N 

W(cu, T, X t ) = v^E - r)) 

OL= 1 

where uj is the test frequency (scale factor) and r is an off¬ 
set relative to the start of the time series (time shift). A 
popular choice for the convolution function (wavelet) is the 
(abbreviated) Morlet wavelet, a plane wave modulated by a 
Gaussian decay profile (determined by the decay constant, 
c): 

f(z) — exp(zcj(t — r) — cuj 2 (t — r) 2 ) 


With unevenly sampled data, however, the DWT is prone 
to spurious high-frequency spikes and fake structures in 
the time-frequency plane. 


Foster 


(1996) proposed that the 
DWT can be reconsidered as a projection of the time series 
onto a trial function (j>(t) = exp— r)) with statisti¬ 
cal weighting w a = exp(—ccc 2 (t — r) 2 ), which compensates 
for irregularities in the data coverage. Furthermore, a bias 
at lower frequencies due to effectively sampling more data 
points through a wider window can be corrected by using Z- 
statistics to give the weighted wavelet Z-transform (WWZ). 

The peak of the WWZ (for constant r) can be used to 
determine the period of any signal in a time series within 
a given time window and the tracks of peaks in the time- 
frequency plane will show any changes of period or ampli¬ 
tude with time. This also allows a check on whether any 
periodicity is present for the whole observed time or is tran¬ 
sient. The resolution of WWZ is dependent on the decay 
constant c used to determine the width of the wavelet win¬ 
dow. If the value is too high then the WWZ power decays 
too rapidly and small amplitude periodic signals will tend 
to be missed whilst if it is too low then the WWZ only has 
broad frequency resolution. A general recommendation is to 
use a value such that the exponential term decreases signif¬ 
icantly in a single cycle Finally, one limitation of the 

WWZ is that it cannot be used for data with gaps larger 
than the period in question. 

An alternate approach finds the wavelets as solutions 
to an eigenvalue problem where they maximize the energy 
in a frequency band \ f — f c \ ^ fw about a central frequency 


Figure 1. The characteristic restframe timescale determined 
from the Slepian wavelet variance as a function of absolute mag¬ 
nitude for the 243,500 quasars in the data set. The black points 
indicate the median value in each magnitude bin. 
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i=/ c . These Slepian wavelets are related to discrete prolate 
spheroidal functions and have the advantage that they can 


be applied to irregular and gappy time series (Mondal & Per- 
cival||20lT ). This type of wavelet analysis makes use of the 


variance of the wavelet coefficients to identify those scales 
making the largest contribution to the total wavelet vari¬ 
ance of a time series. iGraham et al. 20141 used this to iden¬ 
tify a characteristic ~50 day restframe timescale in quasars 
that marks an apparent deviation from the expected behav¬ 
ior for a first-order continuous time autoregressive process 
(CAR(l), also known as a damped random walk or Ornstein- 
Uhlenbeck process). This restframe timescale also strongly 
(anti)correlates with absolute magnitude (see Fig. [TJ . 

Periodic behavior is not expected as a result of a 


CAR(l) process (although see the discussion in Sec. 6.1 


about higher order processes). Periodic objects should there¬ 
fore show a strong deviation in their wavelet variance from 
that expected for a CAR(l) model with a shorter character¬ 
istic timescale than regular quasars. A periodic candidate 
should therefore show as a statistically significant outlier 
from the timescale - absolute magnitude relationship, ap¬ 
pearing below the median curve in Fig. ^ 


2.2 Autocorrelation function (ACF) 

The autocorrelation function (ACF) is a measure of how 
closely a quantity observed at a given time is related to the 
same quantity at another time (the time difference between 
the two is the lag time). In the case of a (quasi-)periodic 
signal, peaks in the ACF will indicate lag times at which 
the signal is perfectly (highly) correlated with itself and 
these can be interpreted as periods at which the signal 
(quasi-)repeats. The peak at zero lag (r = 0) is the strongest 
peak in the ACF of a (quasi-)periodic signal. Successive 
peaks at increasing multiples of the period are increasingly 
weaker, and as the lag value approaches the coherence time 
- the time required for a signal to change its frequency and 
phase information, the peaks in the autocorrelation function 
approach the “noise floor” generated by the non-periodic 
components of the signal and the noise. 


McQuillan et al. (2013) discuss the robustness of the 


ACF for period detection in time series data. Since the ACF 
measures only the degree of self-similarity of the light curve 
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at a given lag time, the period remains detectable even when 
the amplitude and phase of the photometric modulation 
evolve significantly during the timespan of the observations. 
The ACF is also not susceptible to residual instrumental sys- 
tematics since correlated noise, long-term trends and discon¬ 
tinuities manifest as monotonic trends in the ACF, on top of 
which local maxima may still be identified. In fact, the ACF 
is preferable as a period-finding method to the Lomb-Scargle 
periodogram in terms of clarity and robustness. 

For irregularly sampled data, the standard ACF esti¬ 
mator of Edelson & Krolik (1988) can be defined as: 


ACF(kAr) 


Ei=l ti) 

EIU £"=i b k{tj - u) 


where the observations have been normalized to zero mean 
and unit variance before the analysis. The kernel bk(tj — ti) 
selects the observations whose time lag is not further than 
half the bin width from kAr: 


Figure 3. The median absolute deviation - magnitude relation¬ 
ship for the 243,500 quasars in the data set. The black points 
indicate the median value in each magnitude bin. The contours 
indicate the la, 2a and 3a levels. 
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bk(tj - ti) = 1 for | (tj - ti)/Ar - k | < 1/2 

0 otherwise 

However, this estimator has a number of disadvantages in¬ 
cluding a high variance and not always producing positive 
semidefinite covariance matrix estimates (necessary to en¬ 
sure the non-negativity of its Fourier transform; |Rehfeld| 
et al. 2011). Whilst a Gaussian kernel addresses some of 


these, Alexander (2013) recommends Fisher’s ^-transform 


and equal population binning (at least 11 points in each At 
bin) for a more effective estimator. Fig. [2] shows the ACF 
for a sample quasar light curve calculated using both the z- 


transform method and Scargle’s alternative method (Scargle 


1989). 


McQuillan et al. (2013) identify the rotation periods 


of Kepler held M-dwarfs from the largest peaks (relative to 
neighbouring local minima) in their ACFs smoothed with 
a fixed width Gaussian kernel. Our data, however, have 
variable sampling in terms of cadence, temporal baseline 
and number of observations, which makes defining a sin¬ 
gle smoothing width value difficult. We have found that a 
better procedure to determine the period of any signal in 
the data is to fit the ACF of a light curve with a Gaussian 
process model (also known as kriging). Assuming the pro¬ 
cess has a squared-exponential covariance, this provides the 
best linear unbiased prediction of the underlying function 
from noisy data. We define the period as the time-lag value 
associated with the largest peak in the ACF between the sec¬ 
ond and third zero-crossings of the fitted model. This shows 
very good agreement with values determined from ACFs 
smoothed with a range of different width Gaussian kernels. 
We have verified the period determined by this method with 
that obtained with the conditional entropy algorithm of |Gra-| 
ham et al. (2013). We would expect broad agreement be¬ 


tween the period determined from the ACF and that iden¬ 
tified by the wavelet technique. 

Periodically driven stochastic systems are expected to 
have an ACF with the form of an exponentially decaying 
cosine function (Jung 1993). Alternatively, if the light curve 
is the result of a continuous autoregressive moving average 
CARMA process (CARMA) (see Sec. 6.1) then its ACF is 
a sum of exponentially damped sinusoids and exponential 


decays, depending on whether the roots of the characteris¬ 
tic polynomial of the CARMA process are real or imaginary 
(Kelly et al. 2014). To test either case, we have also deter¬ 
mined the best fitting exponentially damped sinusoid to the 
ACF of the form: ACF(r) m Acos(cjt) exp(— At), where the 
amplitude A, period p = 2 tt/uj, and decay A are determined 
using maximum likelihood estimation. 


2.3 Shape and coverage 

In the simplest case, periodicity associated with a Keple- 
rian orbit would show as a pure sinusoidal signal in the 
light curve. Although the projective geometry of a system 
- orbital inclination, eccentricity and argument of perias- 
tron - will distort this, the variation will remain a smooth 
periodic function, essentially a multiharmonic sine function. 
Regular flaring activity and similar phenomena not neces¬ 
sarily associated with a close SMBH binary pair would show 
as (quasi-)periodic discontinuities. We can therefore use the 
scatter around an expected sinusoidal waveform to identify 
those objects most closely exhibiting the desired behavior. 

The best-fit truncated Fourier series (up to 6 terms) is 
determined for a given light curve phased using the period 
calculated by the ACF method (see above): 

6 

<j)(t) — rrio + a n sin(nujt + <p n ) 

n=1 

where (j>(t) is the phased light curve, mo is the median mag¬ 
nitude and uj — 2ir/period acf - An F-test is used to deter¬ 
mine the exact number of harmonics to include. The scatter 
about the best fit will be a combination of the intrinsic vari¬ 
ability of the quasar and the quality of the fit. We therefore 
scale the rms scatter of the phased data about the fit by the 
median absolute deviation (MAD) of the light curve and 
only consider objects with scatter below some cutoff value 
as binary candidates. We note, however, that the MAD is a 
function of magnitude: fainter objects show large MAD val¬ 
ues as there is a larger noise contribution for low S/N (faint) 
sources (see Fig. [3}. We therefore employ a MAD value nor¬ 
malized by the median value for a given magnitude to ensure 
that objects with equivalent variability strength (irrespec¬ 
tive of magnitude) can be compared. 
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Figure 2. The autocorrelation function for the given quasar light curve (left) evaluated using the ^-transform method (middle) and 
Scargle’s method (right). The dashed line shows the best-fit polynomial for the ^-transform fit constrained to go through ACF(O) = 1. 
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Figure 4. The range of SMBH masses and binary separations 
to which a synoptic data set is sensitive. The solid upper line 
for each separation indicates a 2 : = 5 track and the solid lower 
line a z = 0.05 track whilst the two internal dotted lines show 
z = 1.0 (lower) and z = 2.0 (upper) tracks respectively. The 
hatched region indicates the range over which CRTS has temporal 
coverage of 1.5 cycles or more of a periodic signal. The points 
indicate the locations of the CRTS candidates and the best known 
SMBH binary candidate, OJ 287 (solid black star). 





The temporal coverage of the data set used will deter¬ 
mine the range of SMBH masses and binary separations to 
which it is sensitive (see Fig.[4|. A decade-long baseline will 
only sample a fraction of an orbital period for either a low 
mass SMBH pair (where the primary BH mass < 10 7 Mq) 
or a more massive pair at separations > lpc. Since we want 
objects with well-sampled periodic behavior, we limit our 
search to objects with coverage of at least 1.5 cycles, as¬ 
suming their ACF determined period. We also only consider 
objects with more than 50 observations in their light curves 
to minimize any systematic effects in the statistics from poor 
sampling. 


3 DATA SETS 


There are few data sets with sufficient sky and/or tempo¬ 
ral coverage and sampling to support an extensive search 
for quasars exhibiting periodic behavior. Most large stud¬ 
ies of long-term quasar variability, e.g., SDSS with POSS 
( MacLeod et al.||2012 ) or Pan-STARRSl (Morganson et al. 
2014), consist of relatively few epochs of data spread over 


a roughly decadal baseline, which is sufficient to model en¬ 
semble behavior but not to identify specific patterns in in¬ 
dividual objects. CRTS represents the best data currently 
available with which to systematically define sets of quasars 
with particular temporal characteristics. 


3.1 Catalina Real-time Transient Survey (CRTS) 


CRTS leverages the Catalina Sky Survey data streams from 
three telescopes - the 0.7 m Catalina Sky Survey Schmidt 
and 1.5 m Mount Lemmon Survey telescopes in Arizona and 
the 0.5 m Siding Springs Survey Schmidt in Australia - used 
in a search for Near-Earth Objects, operated by Lunar and 
Planetary Laboratory at University of Arizona. CRTS covers 
up to ~2500 deg 2 per night, with 4 exposures per visit, sepa¬ 
rated by 10 min, over 21 nights per lunation. All data are au¬ 
tomatically processed in real-time, and optical transients are 
immediately distributed using a variety of electronic mech- 
anismsj Th e data are broadly calibrated to Johnson V (see 


Drake et al. 2013 


for details) and the full CRTS data sef|^] 
contains time series for approximately 500 million sources. 

The Million Quasars (MQ) catalogu^]v3.7 contains ah 
spectroscopically confirmed type 1 QSOs (309,525), AGN 
(21,728) and BL Lacs (1,573) in the literature up to 2013 
November26 and formed the basis for the results of Graham 
(2015) ^We have extended this with 297,301 spectro- 


et al. 


scopically identified quasars in the SDSS Data Release 12 
( Paris et al.||2015 ). We crossmatched this combined quasar 
list against the CRTS data set with a 3" matching radius 
and find that 334,446 confirmed quasars are covered by the 


3 http://www.skyalert.org 

4 http://nesssi.cacr.caltech.edu/DataRelease 

5 http://quasars.org/milliquas.htm 

6 In this analysis we make no distinction based on object class 
and consider all these sources together. 
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6 M. J. Graham et al. 


Figure 5. The U-band magnitude distribution for the CRTS 
known quasar sample (red). The tail of sources fainter than 
V ~ 20.3 are from the Mount Lemmon 1.5m telescope in the 
Catalina surveys. The stepped distribution (black) shows the re¬ 
jected objects (less than 50 data points) from the sample which 
are primarily fainter sources. 



Magnitude 


Figure 6. The redshift distribution for the CRTS known quasar 
sample (red), omitting the tail of the distribution which contains 
164 quasars with z > 5. The stepped distribution (black) again 
shows the rejected objects. 



full CRTS. Of these, 83,782 do not pass our sampling crite¬ 
rion (i.e., they have less than 50 observations in their light 
curves), leaving a data set of 250,664 quasars. Note that the 
majority of rejected objects are faint {V > 20). The magni¬ 
tude and redshift distributions for our sample are shown in 
Figs. [5] and [6] 


3.2 Preprocessing 

It is common to preprocess data to remove spurious outlier 
points which may be caused by technical or photometric er¬ 
ror. The danger, of course, is removing real signal, although 
a robust method should be unaffected by the presence of 
noisy data. We have created a cleaned version of each data 
set following the procedure of |Palanque-Delabrouille et al.| 
(201lJ (PD11): a 3-point median filter was first applied to 
each time series, followed by a clipping of all points that 
still deviated significantly from a quintic polynomial fit to 
the data. To ensure that not too many points are removed, 
the clipping threshold was initially set to 0.25 mag and then 
iteratively increased (if necessary) until no more than 10% of 
the points were rejected. Note that PD 11 used a 5cr thresh¬ 


old but we found that this removed too few points and so 
used the limit set by Schmidt et al. (2010) instead. 

There is always a chance that the light curve for a 
particular object may be contaminated by stray light (e.g., 
diffraction spikes, reflections) from nearby bright sources (up 
to ~2' away for CRTS), low surface brightness galaxies or 
genuine blends, which can give a false seeing-dependent re¬ 
sult in the analysis. If an object has a nearby counterpart 
(within ~5 /x ) which is unresolved by CRTS (~2 pixels), the 
light curve will be the combined signal from both, dominated 
by the brighter of the two. We have defined a distance and 
magnitude-based criterion to exclude quasars near bright 
stars by crossmatching the quasar list against the Tycho cat¬ 
alog of bright stars (V < 14) and reject 5,319 sources. We re¬ 
ject all objects within a 350" radius of any star brighter than 
V = 6.67 and within a radius of 200 // /(R — 6) for V > 6.67. 
To identify blends, we crossmatched our list against SDSS 
and flagged those objects where there is a significant dif¬ 
ference (taking into account the variability of the source) 
between the magnitude of the CRTS source and the SDSS r 
of its match, typically indicating a faint close SDSS compan¬ 
ion contributing to the CRTS photometry: \V — r\ > 5 &v- 
We also flagged those CRTS quasars outside the SDSS foot¬ 
print (or without a SDSS companion) where the distribution 
of the spatial positions of individual photometry points in 
the light curve show significant scatter, e.g., a bimodal dis¬ 
tribution indicating two sources: a pos > 5 a p os(V) where 
a p os(V) is the expected positional scatter for all quasar of 
magnitude V. Together these reject another 4.560 sources. 
The final size of the data set is 243.486 quasars. 


3.3 Mock data 


To first approximation, quasar variability is well-described 
by a Gaussian first-order continuous autoregressive model 
(CAR(l)), also known as a damped random walk or 
Ornstein-Uhlenbeck process (see |Graham et al. 2014| [Kelly 
|et al. 20l4| for details of where this description breaks down). 
Formally, the temporal behavior of the quasar flux X(t) is 
given by: 

dX(t) — -—X(t)dt-\- aVdte(t) + bdt , r,a,t> 0 

T 


where r is the relaxation time of the process, a is the vari¬ 
ability of the time series on timescales short compared to r, 
br is the mean magnitude, and e(t) is a white noise process 
with zero mean and variance equal to 1. For each quasar 
in our data set, we have generated a simulated light curve 
assuming that it follows a CAR(l) model. Using the actual 
observation times, U, we replace the observed magnitudes 
with those that would be expected under a CAR(l) model. 
The magnitude X(t) at a given timestep At from a previous 
value X(t — At) is drawn from a Gaussian distribution with 
mean and variance given by, e.g., ( Kelly et al.|2009 ): 


E(X(t)\X(t - At)) = e~ At,T X(t - At) + 6r(l - e~ At/r ) 


V & v(X(t)\X(t - At)) = T X {1 _ e - 2 **/r] 


We add a Gaussian deviate normalized by the photometric 
error associated with the magnitude to be replaced at each 
time t to incorporate measurement uncertainties into the 
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mock light curves. For each light curve, we set br to its me¬ 
dian value and use the rest frame CAR(l) fitting functions 
determined by MacLeod et al. | ((2010): 


los/ = +Bloe I sss) + c(Jt+23) + D '°* (mS~ e 


where (A,B,C,D) = (-0.51,-0.479,0.131,0.18) for / = 
SFov and (A,B,C,D) = (2.4,0.17,0.03,0.21) for / = r. 
Note that a value of C — 0.113 is used for / = SF a 0 for 
non K-corrected magnitudes. Mi is the absolute magnitude 
of the quasar and A rf is the restframe wavelength of the 
filter - here taken to be SDSS r, A rf — 6250A/(1 + z). 
The mass of the black hole given the absolute magnitude is 
drawn from a Gaussian distribution: 


p( log Mbh\Mi) 


1 

V / 27tct 


exp 


(log M B h - p ) 2 
2cr 2 


where p = 2.0 — 0.27 Mi and a — 0.58 + 0.011 Mi. This is 
based on Shen et al. (2008) results. Differences in cosmolo¬ 
gies used in estimating the best-fit parameters for CAR(l) 
and black hole mass should only have a 1% effect (MacLeod 
thesis, 2012). 


4 RESULTS 

We have applied our technique which identifies strong peri¬ 
odic behavior to both the real CRTS quasar data set and its 
simulated counterpart. From the real data set, we define an 
initial candidate sample according to the following criteria 
(see Table [I] for a summary). For the wavelet component, we 
set the decay constant c — 0.001 and use the frequency range 
0.0003 - 0.05 day -1 . As we are interested in well-defined pe¬ 
riodicity, we only consider those objects whose wavelet peak 
significance places them in the top quartile of the data set 
(in terms of significance), which translates to a WWZ peak 
value cutoff of wwzpk > 50. We also only consider objects 
with a characteristic timescale from their Slepian wavelet 
variance at least la below the expected value for their ab¬ 
solute magnitude. This reflects that periodic objects are not 
expected to be well-described by a CAR(l) model. 

The first ACF-based constraint is that the periods de¬ 
termined by the WWZ and the ACF, respectively, should 
agree to within 10 per cent. The shape of the ACF should 
also be reasonably described by an exponentially damped 
cosine with not too small an amplitude (A > 0.3) and a de¬ 
cay constant such that the ACF drops by at most a factor 
of 1/e over the temporal baseline of the time series. This 
corresponds to A < 10 -3 for the shortest time coverage in 
the data set and so we use this as a fiducial value. 

In terms of the shape constraint, we limit the selection 
to those objects whose rms scatter about the best-fit trun¬ 
cated Fourier series to their phased light curve is less than 
the la lower limit on the magnitude-normalized median ab¬ 
solute deviation, i.e., rms/mad < 0.67. Finally, as previ¬ 
ously described, we also limit the allowed temporal coverage 
to 1.5 cycles or greater, i.e. t/ period acf > 1.5, where r is 
the timespan of a light curve, and only consider light curves 
with 50 or more observations. 

Ideally, we want to determine the discriminating hyper¬ 
plane in the multidimensional feature space that provides a 
clear separation between periodic and non-periodic objects. 


The feature cuts that we have defined provide a first approxi¬ 
mation to this but likely candidates close to a selection value 
may have been missed due to error, poor sampling, etc. We 
therefore use the initial sample defined by the feature cuts as 
part of the training set for a support vector machine (SVM) 
with a radial basis function kernel. A number of similar ma¬ 
chine learning algorithms could also be employed here but 
we found from simulations that the SVM was the best in 
terms of speed and accuracy. 

The other part of the training set comprises a set of 
randomly selected quasars that do not meet the periodic¬ 
ity selection criteria. From simulations with known num¬ 
bers of periodic objects and periodic/non-periodic ratios, we 
have determined that the most adequate training set mix is 
~100:1, given our expected detection rate (see Section 6.3). 
10000 randomly selected “non-periodic” quasars are thus 
added to the training set. The trained SVM is then used 
to classify each quasar in the full data set as periodic or 
not periodic based on their feature values. We repeat this 
process 50 times and identify those objects in the full data 
set that are selected each time as periodic as a final binary 
candidate. 

Our final sample consists of 111 quasars (see Table [5] for 
details and Fig. |]for their light curves). We have also applied 
the same selection process to the simulated data set of ob¬ 
jects with CAR(l) model light curves and find that none of 
them satisfy our criteria. If we apply the SVM trained on the 
real data to the simulated data set, we also find that that no 
objects are selected in ah iterations as periodic. There are, 
however, three mock sources which are selected in almost 
ah iterations (see Fig.[7|. Though these objects are the clos¬ 
est to our selection criteria, they look qualitatively different 
from the real sample with no clear periodicity. The range of 
variability they exhibit is also significantly more than in the 
real sample: their mock light curves have MAD values two 
or three times the values seen in the real data. 

If we take the simulated data results as an indication of 
the expected number of quasars to show comparable periodic 
behavior by chance then the number of real objects selected 
is statistically significant, particularly with the higher tem¬ 
poral coverage constraint. It also suggests strong periodicity 
(on an approximately decadal timescale) is not common in 
quasars nor is it expected as an artifact of a CAR(l) pro¬ 
cess. We note, however, that our selection criteria will select 
against longer periodic as well as more general quasi-periodic 
behavior as may be expected from SMBH binaries at greater 
than 1 pc separation. 


5 BINARY CANDIDATES IN THE 
LITERATURE 

5.1 Spectroscopic 

As previously mentioned, the availability of multiepoch 
spectroscopy enables searches for SMBH binaries on the ba¬ 
sis of time-dependent velocity shifts in the broad lines, either 
with respect to each other or narrow lines associated with 
the host galaxy. It is interesting to see whether any spectro¬ 
scopically selected candidates in the literature show periodic 
photometric variability in CRTS data. We have examined 
the light curves of 30 spectroscopic SMBH binary candi¬ 
dates from the literature (see Table [3| . We only consider 
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Table 1. The selection constraints used to define an initial sample of binary candidates and the respective counts for the real and 
simulated data sets for each constraint. 


Component 

Constraint 

Real 

Mock 

Wavelet peak value 

wwzpk > 50 

24437 

32078 

Slepian wavelet deviation 

Tslep < t(M v ) - cr T 

37828 

27746 

WWZ-ACF period 

0.9 < Pacf/pwwz < 1-1 

30330 

63637 

ACF amplitude 

A > 0.3 

108625 

143447 

ACF decay 

A < 10" 3 

172049 

182592 

Shape 

rms/MAD < 0.67 

11794 

3944 

Temporal coverage 

r/pACF > 1.5 

245234 

182257 

Number of points 

n > 50 

243486 

243486 

Combined 

_ 

101 

0 

Final 

- 

111 

0 


Figure 7. The three objects in the simulated data which most closely match the selection criterion. The best-fit sinusoid for the 
determined period for each source is plotted in red. 





here the most likely binary candidates reported, e.g., classi¬ 


fied as such from the shape of their Balmer lines (Tsalmantza 


2011 Decarli et al. 2013), and not secondary candidates also 


reported with asymmetric line profiles, double-peaked emit¬ 
ters, or other features since these do not sufficiently indicate 
a SMBH binary. Although the predicted separation of these 
candidates is greater than in our candidates, it is still typi¬ 
cally less than 1 pc and we would expect there to be some 
detectable behavior in the photometric data. 

There seem to be a number of morphologically distinct 
behaviors seen with these candidates. To quantify this, we 
have defined a number of template functions that capture 
different shape profiles: 

Linear, a line with gradient given by the Theil-Sen estimator 
for the light curve which passes through the median data 


point {trned’i'hTlmed) 5 

Polynomial : the best fitting polynomial up to 10 th order 
determined by a F-test between successive orders; 
Sinusoidal : the best fitting sinusoid with period determined 
by the ACF method (see section 2.2); 

Logistic: a generalized logistic function fit to the light curve: 


m(t) = A + 


K -A 

(1 + 


Hump/Dip: a linear least-squares fit to the light curve plus 
a quadratic to the largest contiguous amount of data above 
or below the linear fit. 

For each light curve of interest, we have determined 


which template function provides the best chi-squared fit. 
The shape of the light curve may also just be a random pat¬ 
tern arising from the stochastic nature of the variability and 
so we have also determined both the best-fit CAR(l) and 
CARMA(p,q) model parameters for each light curve and 
checked for correlations between these and the best-fitting 
template. We note that none of the spectroscopic candidates 
shows the type of strong periodicity exhibited by the candi¬ 
dates identified in this paper. 


The division of template shapes is fairly uniform, al¬ 
though objects from Liu et al. (2014) tend to be best-fit 
by the hump template. These candidates specifically show 
significant radial accelerations in their broad H (3 lines from 
dual epoch spectroscopy whereas other studies are based on 
more or other lines. There is also some suggestion of prefer¬ 
ential localization in the CAR(l) r — a 2 plane with objects 
best fit by a sinusoid function lying outside the la contour. 
Those best fit by CARMA models with p — 7 also show 
the same effect but to a higher degree. However, the sample 
size here is too small to determine whether these effects are 
statistically significant but they certainly warrant further 
study. 

Liu, Shuo h Komossa (2014) proposed SDSS 


J120136.02+300305.5 as the first milliparsec SMBH binary 
in a quiescent galaxy on the basis of a candidate tidal dis¬ 
ruption event in its X-ray light curve. Its optical light curve 
shows no significant variation - it has a MAD value of 0.07 
mag which is the median value for its magnitude - or evi- 
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Table 2. The 111 periodic quasar candidates meeting the selection criteria. The full version is available in electronic form online. Black 
hole masses are obtained from Shen et al. (2008) and updates or directly from spectra where a value is not available. Binary system 
parameters - separation r, the restframe merger time ( ti nsp ) and the maximum induced timing residual amplitude (At gw) ~ all assume 
q = 0.5. 


Id 

RA 


Dec 



z 

^median 

Period 

(days) 


r 

(pc) 

tinsp 

(yrs) 


Ataw 

(ns) 

UM 211 

00 

12 

10.9 

-01 

22 

07.6 

1.998 

17.38 

1886 

- 

- 

- 



- 

UM 234 

00 

23 

03.2 

+01 

15 

33.9 

0.729 

17.82 

1818 

9.19 

0.013 

1.2 

X 

10 4 

3.0 

SDSS J014350.13+141453.0 

01 

43 

50.0 

+14 

14 

54.9 

1.438 

17.68 

1538 

9.21 

0.009 

2.8 

X 

10 3 

2.4 

PKS 0157+011 

02 

00 

03.9 

+01 

25 

12.6 

1.170 

18.04 

1052 

- 

- 

- 



- 

RX J024252.3-232633 

02 

42 

51.9 

-23 

26 

34.0 

0.680 

19.01 

1818 

- 

- 

- 



- 

US 3204 

02 

49 

28.9 

+01 

09 

25.0 

0.954 

18.09 

1666 

8.95 1 

0.009 

1.7 

X 

10 4 

1.0 

CT 638 

03 

18 

06.5 

-34 

26 

37.4 

0.265 

16.68 

1515 

- 

- 

- 



- 

RXS J04117+1324 

04 

11 

46.9 

+13 

24 

16.5 

0.277 

16.20 

1851 

8.16 

0.006 

1.4 

X 

10® 

0.1 

HS 0423+0658 

04 

26 

30.2 

+07 

05 

30.3 

0.170 

15.73 

1123 

- 

- 

- 



- 

2MASS J04352649-1643460 

04 

35 

26.5 

-16 

43 

45.7 

0.098 

16.96 

1369 

7.78 

0.004 

4.1 

X 

10® 

0.1 

SDSS J072908.71+400836.6 

07 

29 

08.6 

+40 

08 

37.0 

0.074 

16.08 

1612 

5.71 

0.001 

1.9 

X 

10 10 

0.0 

SDSS J080237.60+340446.3 

08 

02 

37.6 

+34 

04 

46.6 

1.119 

18.40 

1428 

8.96 

0.007 

8.8 

X 

10 3 

0.9 

SDSS J080648.65+184037.0 

08 

06 

48.6 

+18 

40 

37.3 

0.745 

20.61 

0892 

7.99 

0.003 

1.7 

X 

10 5 

0.0 

SDSS J080809.56+311519.1 

08 

08 

09.5 

+31 

15 

18.9 

2.642 

18.91 

1162 

8.36 

0.003 

1.2 

X 

10 4 

0.1 

SDSS J081133.43+065558.1 

08 

11 

33.4 

+06 

55 

58.3 

1.266 

18.48 

1587 

9.39 

0.010 

1.8 

X 

10 3 

5.0 

SDSS J081617.73+293639.6 

08 

16 

17.8 

+29 

36 

40.7 

0.768 

17.80 

1162 

9.77 

0.013 

3.6 

X 

10 2 

23.7 

FBQS J081740.1+232731 

08 

17 

40.2 

+23 

27 

32.0 

0.891 

16.58 

1190 

9.55 

0.011 

7.4 

X 

10 2 

9.6 

SDSS J082121.88+250817.5 

08 

21 

22.0 

+25 

08 

16.2 

1.906 

17.74 

1886 

9.53 

0.011 

8.8 

X 

10 2 

8.2 

SDSS J082716.85+490534.0 

08 

27 

16.9 

+49 

05 

34.9 

0.682 

18.38 

1612 

8.96 

0.009 

2.2 

X 

10 4 

1.3 

SDSS J082827.84+400333.9 

08 

28 

27.8 

+40 

03 

34.1 

0.968 

17.79 

1886 

8.87 

0.008 

3.1 

X 

10 4 

0.8 

SDSS J082926.01+180020.7 

08 

29 

26.0 

+18 

00 

20.7 

0.810 

19.81 

1449 

8.42 

0.005 

1.1 

X 

10 5 

0.1 

SDSS J083349.55+232809.0 

08 

33 

49.6 

+23 

28 

09.2 

1.155 

17.58 

1086 

9.40 

0.008 

7.4 

X 

10 2 

4.7 

SDSS J084146.19+503601.1 

08 

41 

46.3 

+50 

36 

00.5 

0.555 

18.49 

1694 

7.44 

0.003 

1.1 

X 

10 7 

0.0 

BZQJ0842+4525 

08 

42 

15.3 

+45 

25 

45.0 

1.408 

17.20 

1886 

9.48 

0.012 

1.7 

X 

10 3 

7.3 

SDSS J091554.50+352949.6 

09 

15 

54.5 

+35 

29 

49.9 

0.896 

17.92 

1369 

9.05 

0.008 

7.3 

X 

10 3 

1.5 

SBS 0920+590 

09 

23 

58.7 

+58 

49 

06.3 

0.709 

16.76 

0649 

8.76 

0.004 

4.1 

X 

10 3 

0.4 

SDSS J092911.35+203708.5 

09 

29 

11.3 

+20 

37 

09.2 

1.845 

18.51 

1785 

9.92 

0.014 

1.8 

X 

10 2 

36.2 

HS 0926+3608 

09 

29 

52.1 

+35 

54 

49.6 

2.150 

16.99 

1562 

9.95 

0.012 

8.4 

X 

10 1 

38.8 

SDSS J093819.25+361858.7 

09 

38 

19.3 

+36 

18 

58.9 

1.677 

18.80 

1265 

9.32 

0.007 

8.3 

X 

10 2 

3.3 

SDSS J094450.76+151236.9 

09 

44 

50.7 

+15 

12 

37.5 

2.118 

17.72 

1428 

9.61 

0.009 

2.6 

X 

10 2 

10.0 

SDSS J094715.56+631716.4 

09 

47 

15.6 

+63 

17 

17.3 

0.487 

16.10 

1724 

9.22 

0.012 

1.3 

X 

10 4 

4.3 

HS 0946+4845 

09 

50 

00.7 

+48 

31 

29.9 

0.589 

16.87 

1587 

8.59 

0.007 

1.0 

X 

10 5 

0.3 

KUV 09484+3557 

09 

51 

23.9 

+35 

42 

49.2 

0.398 

17.77 

1162 

8.31 

0.005 

1.8 

X 

10 5 

0.1 

SDSS J102255.21+172155.7 

10 

22 

55.2 

+17 

21 

56.0 

1.062 

18.48 

1666 

8.69 

0.006 

3.9 

X 

10 4 

0.4 

SDSS J102349.38+522151.2 

10 

23 

49.5 

+52 

21 

51.8 

0.955 

16.96 

1785 

9.59 

0.014 

1.8 

X 

10 3 

12.1 

RXS J10304+5516 

10 

30 

25.0 

+55 

16 

23.4 

0.435 

16.74 

1515 

8.43 

0.006 

2.2 

X 

10 5 

0.2 

SDSS J103111.52+491926.5 

10 

31 

11.5 

+49 

19 

27.2 

1.203 

17.64 

1612 

9.04 

0.008 

7.8 

X 

10 3 

1.3 

SDSS J104430.25+051857.2 

10 

44 

30.3 

+05 

18 

56.8 

0.905 

17.33 

1333 

9.24 

0.009 

3.3 

X 

10 3 

3.0 

SDSS J 104758.34+284555.8 

10 

47 

58.3 

+28 

45 

56.2 

0.792 

19.09 

1851 

8.72 

0.008 

6.9 

X 

10 4 

0.5 

SDSS J104941.01+085548.4 

10 

49 

41.0 

+08 

55 

48.5 

1.185 

17.87 

1428 

9.37 

0.009 

1.7 

X 

10 3 

4.5 

MS 10548-0335 

10 

57 

22.3 

-3 51 31.3 

0.555 

17.67 

0892 

- 

- 

- 



- 

CSO 67 

11 

03 

27.5 

+29 

48 

11.2 

0.909 

17.52 

2083 

9.35 

0.013 

7.1 

X 

10 3 

5.2 

SDSS J 110554.78+322953.7 

11 

05 

54.8 

+32 

29 

54.1 

0.151 

17.45 

1724 

8.24 

0.007 

1.2 

X 

10® 

0.2 

SDSS J113050.21+261211.4 

11 

30 

50.2 

+26 

12 

11.8 

1.012 

17.69 

2173 

9.32 

0.013 

7.6 

X 

10 3 

4.6 

SDSS J 113916.47+254412.6 

11 

39 

16.4 

+25 

44 

13.0 

1.012 

17.61 

2439 

9.16 

0.012 

1.9 

X 

10 4 

2.6 

SDSS J114438.34+262609.4 

11 

44 

38.3 

+26 

26 

10.1 

0.974 

18.39 

1315 

9.38 

0.010 

1.7 

X 

10 3 

4.9 

SDSS J114749.70+163106.7 

11 

47 

49.7 

+16 

31 

06.8 

0.554 

18.78 

1449 

8.22 

0.005 

3.5 

X 

10 5 

0.1 

SDSS J114857.33+160023.1 

11 

48 

57.4 

+16 

00 

22.7 

1.224 

18.14 

1851 

9.90 

0.017 

4.1 

X 

10 2 

37.6 

SDSS J115141.81+142156.6 

11 

51 

41.8 

+14 

21 

57.0 

1.002 

18.11 

1492 

9.11 

0.009 

6.3 

X 

10 3 

1.8 

SDSS J115346.39+241829.4 

11 

53 

46.4 

+24 

18 

29.8 

0.753 

18.41 

1666 

8.97 

0.009 

2.1 

X 

10 4 

1.3 

SDSS J 121018.34+015405.9 

12 

10 

18.3 

+01 

54 

06.2 

0.216 

16.90 

1612 

8.54 

0.008 

2.6 

X 

10 5 

0.6 

SDSS J121018.66+185726.0 

12 

10 

18.7 

+18 

57 

27.0 

1.516 

17.42 

1754 

9.53 

0.011 

1.1 

X 

10 3 

8.4 

SDSS J 121056.83+231912.5 

12 

10 

56.8 

+23 

19 

13.0 

1.260 

18.47 

1785 

8.78 

0.007 

2.7 

X 

10 4 

0.5 
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Table 2 - continued 


Id 

RA 


Dec 

z 

^median 

Period 

(days) 


r 

(pc) 

tinsp 

(yrs) 


Atcw 

(ns) 

SDSS J121457.39+132024.3 

12 

14 

57.4 

+13 20 24.5 

1.494 

18.59 

1923 

9.46 

0.011 

1.8 

X 

10 3 

6.6 

SDSS J123147.27+101705.3 

12 

31 

47.3 

+10 17 05.4 

1.733 

18.83 

1851 

9.20 

0.009 

3.5 

X 

10 3 

2.3 

SDSS J123821.84+030024.2 

12 

38 

21.8 

+03 00 24.6 

0.380 

18.46 

1250 

8.92 

0.008 

2.2 

X 

10 4 

1.4 

SDSS J 124044.49+231045.8 

12 

40 

44.5 

+23 10 46.1 

0.722 

18.40 

1428 

8.94 

0.008 

1.6 

X 
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Figure 8. Light curves for periodic candidates. CRTS data (black) is shown; complementary data (blue) from the LINEAR survey ( |Sesar| 
|et a l. 2011) is also included where available. The best-fit sinusoid for the determined period is plotted in red. The full set is available 
online. 
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Figure 8 - continued 
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Figure 8 - continued 
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Figure 8 - continued 
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Figure 8 - continued 
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Figure 8 - continued 
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Figure 8 - continued 
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Figure 8 - continued 
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Figure 8 - continued 
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Figure 8 - continued 
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Table 3. SMBH binary candidates identified in the literature 


Id 

RA 

Dec 

Shape 

Best fits 
Poly log 10 r 

logio cr2 

CARMA 

MAD 

Source 

J0049+0026 

00 49 19.0 

+00 26 09.4 

Hump 

3 

2.999 

-4.579 

(2, 0) 

0.10 

Optical 10 

J0301+0004 

03 01 00.2 

+00 04 29.3 

Periodic 

1 

-1.957 

0.753 

(4, 0) 

0.21 

Optical 9 

J0322+0055 

03 22 13.9 

+00 55 13.4 

Polynomial 

5 

2.720 

-4.327 

(6, 0) 

0.07 

Optical 9 

J0757+4248 

07 57 00.7 

+42 48 14.5 

Logistic 

3 

3.235 

-5.031 

(2, 0) 

0.11 

Optical 10 

J0829+2728 

08 29 30.6 

+27 28 22.7 

Polynomial 

4 

3.095 

-3.827 

(4, 0) 

0.23 

Optical 8,12 

J0847+3732 

08 47 16.0 

+37 32 18.1 

Polynomial 

4 

2.844 

-4.001 

(4, 0) 

0.12 

Optical 8 

J0852+2004 

08 52 37.0 

+20 04 11.0 

Hump 

4 

2.594 

-3.871 

(4, 0) 

0.15 

Optical 8 

OJ 287 

08 54 48.9 

+20 06 30.6 

Periodic 

5 

1.678 

-2.021 

(7, 0) 

0.34 

Optical 1 

J0927+2943 

09 27 04.4 

+29 44 01.3 

Periodic 

2 

2.882 

-4.639 

(1, o) 

0.09 

Optical 2,7,11 

J0928+6025 

09 28 38.0 

+60 25 21.0 

Periodic 

1 

2.650 

-4.105 

(6, 0) 

0.10 

Optical 8 

J0932+0318 

09 32 01.6 

+03 18 58.7 

Polynomial 

5 

2.921 

-4.119 

(4, 0) 

0.13 

Optical 6,7,11 

J0935+4331 

09 35 02.5 

+43 31 10.7 

Hump 

5 

2.946 

-4.541 

(6, 0) 

0.05 

Optical 10 

4C 40.24 

09 48 55.3 

+40 39 44.6 

Polynomial 

5 

2.900 

-4.500 

(2, 0) 

0.09 

Radio 13 

J0956+5350 

09 56 56.4 

+53 50 23.2 
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1 

2.713 

-4.392 

(4, 0) 

0.10 
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J1000+2233 

10 00 21.8 

+22 33 18.6 

Polynomial 
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2.142 

-3.593 

(7, 0) 

0.12 

Optical 5,7,11 

J1012+2613 

10 12 26.9 

+26 13 27.2 

Logistic 

2 

2.936 

-4.557 

(5, 0) 

0.09 

Optical 7,11 

J1030+3102 
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3 
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-4.134 
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J1050+3456 
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Polynomial 

10 

2.689 

-4.154 

(2, 0) 

0.12 

Optical 4,7,11 

J1100+1709 

11 00 51.0 

+17 09 34.3 

Hump 

1 

2.906 

-4.232 

(6, 0) 

0.12 

Optical 8,12 

J1112+1813 

11 12 30.9 

+18 31 11.4 

Polynomial 

6 

2.784 

-4.264 

(3, 0) 

0.09 

Optical 8 

J1154+0134 

11 54 49.4 

+01 34 43.6 

Logistic 

1 

2.906 

-3.853 

(5, 2) 

0.15 

Optical 7,11 

J1229-0035 

12 29 09.5 

-00 35 30.0 

Polynomial 

4 

2.834 

-3.947 

(6, 5) 

0.13 

Optical 9 

J1229+0203 

12 29 06.7 

+02 03 08.7 

Polynomial 

5 

1.397 

-3.057 

(7, 2) 

0.07 

Radio 13 

J1305+1819 

13 05 34.5 

+18 19 32.9 

Polynomial 

8 

2.728 

-4.899 

(5, 0) 

0.04 

Optical 8,12 

J1345+1144 

13 45 48.5 

+11 44 43.5 

Hump 

4 

2.440 

-4.039 

(6, 0) 

0.08 

Optical 8 

J1410+3643 

14 10 20.6 

+36 43 22.7 

Logistic 

3 

3.240 

-4.434 

(5, 0) 

0.15 

Optical 9 

J1536+0441 

15 36 36.2 

+04 41 27.1 

Polynomial 

10 

2.732 

-4.474 

(6, 0) 

0.07 

Optical 3,7,11 

J1537+0055 

15 37 06.0 

+00 55 22.8 

Polynomial 

2 

3.105 

-5.987 

(2, 0) 

0.03 

Optical 9 

J1539+3333 

15 39 08.1 

+33 33 27.6 

Periodic 

1 

3.042 

-5.591 

(2, 0) 

0.06 

Optical 7,11 

J1550+0521 

15 50 53.2 

+05 21 12.1 

Hump 

1 

2.804 

-4.817 

(2, 0) 

0.03 

Optical 9 

J1616+4341 

16 16 09.5 

+43 41 46.8 

Logistic 

1 

-1.383 

-0.033 

(4, 0) 

0.18 

Optical 10 

J1714+3327 

17 14 48.5 

+33 27 38.3 

Polynomial 

7 

3.060 

-4.865 

(1, o) 

0.06 

Optical 7,11 

PKS 2155-304 

21 58 52.1 

-30 13 32.1 

Logistic 

1 

2.483 

-2.648 

(3, 1) 

0.44 

Radio 13 

3C 454.3 

22 53 57.7 

+16 08 53.6 

Polynomial 

4 

1.978 

-1.817 

(6, 0) 

0.55 

Radio 13 

1ES 2321+419 

23 23 52.1 

+42 10 59 

Polynomial 

5 

2.531 

-2.974 

(6, 0) 

0.29 

X-ray 14 

J2349-0036 

23 49 32.8 

-00 36 45.8 

Periodic 

1 

-0.337 

-1.153 

(7, 3) 

0.10 

Optical 9 


References: [1] Valtonen et al. (2008); [2] Dotti et al. (2009); [3] Boroson & Lauer 2009; [4] Shields, Bonning & Salviander (2009) ; [5] 
Decarli et al. (2010); [6] Barrows et al. (2011); [7] Tsalmantza et al. (2011); [8] Liu et al. (2014); [9] Shen et al. (2013); [10] |Ju et al.| 
(2 013] ); [11] |Decarli et a l. 2013] [12] Eracleous et al. (2012); [13] Ciaramella et al. (2004); [14] Rani et al. (2009) 


dence of an associated flare and there is no structure to it. 
This suggests that this population of SMBH binaries will 
only be detected in the optical as a result of serendipity 
(assuming that the candidate is an actual binary). 


5.2 Blazars 


Similarly we have examined the CRTS light curves of 
the few blazars reported in the literature as exhibiting 
(quasi-)periodic behavior, e.g., PKS 2155-304 and 1ES 
2321+419 (Rani, Wiita & Gupta]2009). We note that these 
do not show the consistent periodicity that we see in the 
candidates in this analysis. However, their CAR(l) and 
CARMA fits are comparable with those seen for optical can¬ 
didates in the literature as in the previous section. 


5.3 Photometric 


Liu et al. (2015) report the discovery of a periodic quasar 


candidate, FBQS J221648.7+012427, out of a sample of 
168 quasars in the Pan-STARRSl Medium Deep Survey 
(PSMDS) with an average of 350 detections in four filters. It 
has an observed period of 542 days, a redshift of 2.06, and an 
inferred separation of 0.006 pc. Whilst this object is part of 
the data we have considered, it does not appear as a candi¬ 
date in our analysis. Assuming the quoted period, its CRTS 
light curve covers 5.8 cycles but no significant or consistent 
periodic signal is identified by either the WWZ or ACF al¬ 
gorithms and its Slepian wavelet characteristic timescale is 
not a statistically significant outlier. 

The candidate was identified using a Lomb-Scargle (LS) 
periodogram to look for evidence of periodicity. A similar 
approach was used by MacLeod et al. (2010| in a search of 
8863 SDSS Stripe 82 quasar light curves (with an average 
of over 60 observations) with 88 candidates identified. None 
were considered significant, however, as the standard likeli- 
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Figure 8 - continued hood calculation assumes white noise as the null hypothesis 

rather than colored noise as appropriate for an underlying 
CAR(l) process. The technique is regarded as not very pow¬ 
erful for poorly sampled light curves. It is also interesting to 
note that both analyses have a detection rate of about 1 in 
100 rather than the far more conservative 1 in 10000 that we 
find, which is much more in line with expected rates from 
population models (see section [R3| . 

The PSMDS candidate FBQS J221648.7+012427 is also 
reported to have a (restframe) inspiral time of ~ 7 years 
thus putting its final coalescence within the timeframe of 
pulsar timing array experiments within the next decade or 
0 500 1000 1500 2000 2500 3000 350©o. Given that the merger timescale is typically 10 7 years, it 

MJD - 53500 seems statistically unlikely that an object would be found in 

its final phase in such a small sample size. The inspiral time 
is, however, very dependent on the black hole mass and so 
it is possible that the authors have overestimated the mass 
of the system, particularly as the line and continuum fluxes 
were not determined from an electronic format spectrum. 


SDS5 J224829.47 + 1444L8.0 



6 DISCUSSION 


6.1 Quasiperiodic behavior in stochastic models 

Although CAR(l) models have been used as popular statis¬ 
tical descriptors of quasar optical variability, e.g., Ke lly et| 
al. ( 2009| ; [MacLeod et al. (2012), there is a growing body 
of evidence that a more sophisticated approach is required 


( Mushotzky et al.|2011 Zu et al.|2013 Graham et al. 2014 


Kasliwal, Vogeley & Richards 2015). Kelly et al. 2014 de¬ 


scribe the use of CARMA models as a better method to 
quantify stochastic variability in a light curve. (Note that a 
CAR(l) model is the same as a CARMA(1, 0) model). A 
zero-mean CARMA(p,q) process y{t) is defined to be the 
solution to the stochastic differential equation: 


d?y(t) d p ~ 1 y(t) 

h a p -1 — 7| _,-h ■ 


dtp 


A 


dtP~ 


+ a 0 y(t) = 


d q e{t) d q 1 e{t) 

I- P 9-1 j.„-i + • • • + e(t) 


dt q dt q ~ 

where e(t) is a continuous time white noise process with zero 
mean and variance a 2 . This is associated with a character¬ 
istic polynomial: 


^) = E 


a k z 


with roots n,..., r p . The power spectral density (PSD) for 
a CARMA(p, q) process is given by: 

H I I 2 

Now A(z) can be written in terms of its roots: 


a o= n c-rR =o 

3 = 1 

where r = a + ib and so the PSD is proportional to: 
P(w) oc 1 


nu d + iw- bj ) 2 
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For complex roots (and p > 1), local maxima will exist 
at w = bj which correspond to quasi-periodicities in the 
data. If the roots have negative real parts and q < p then 
the CARMA process is stationary. By analogy with under¬ 
damped second order systems, the quasi-periodicities can be 


characterized by the damping ratio, Q = —a,j/ a 2 + b 2 : for 
a given frequency, bj . Smaller values of £ indicate a stronger 
oscillating signal (less damping) with £ = 0 defining a pure 
sinusoidal signal. The bandwidth of the signal, i.e., the range 
of frequencies associated with it, is given by A/ = £5/7r. 

For each of our periodic candidates, we have determined 
the best- fitting CA RM A(p,q ) model using the carma_pack 


cod^] of Kelly et al. 2014 


and checked whether there is 
any quasi-periodic behavior expected at its identified period 
from peaks in the PSD. Out of the 111 candidates, we find 
that five show a detected periodicity which overlaps with 
that expected from its best-fit CARMA model. However, 
the bandwidth of all of these is such that any period within 
the temporal baseline of the object would fit. None of the pe¬ 
riodic behavior is therefore consistent with that potentially 
expected from stochastic variability. 

We have also generated 100 mock light curves of the 
appropriate CARMA(p,q) type for each periodic candidate 
with the same time sampling as the CRTS data and deter¬ 
mined the expected wavelet and ACF parameters. The me¬ 
dian differences between the real value and CARMA model- 
based values are statistically significant in all cases. From 
this we conclude that CARMA(p,q) models do not provide 
an adequate statistical description for our periodic candi¬ 
dates as they cannot reproduce the same range of features 
values seen. 


6.2 Physical interpretation 

Although we have searched for optical periodicity in quasars 
on the assumption that this is associated with a close SMBH 
binary system, the exact physical mechanism responsible for 
this behavior is uncertain. The sinusoidal nature of the sig¬ 
nal suggests that it is kinematic in origin and there are a 
number of potential explanations that should be considered. 


6.2.1 A jet with a varying viewing angle 


The optical flux could be the superposition of thermal emis¬ 
sion from the accretion disk and a differential Doppler 
boosted non-thermal contribution associated with a jet 


(Rieger 2004). The observed periodicity would be the re¬ 


sult of a periodically varying viewing angle, driven either 
by the orbital motion of a SMBH binary system, a precess- 
ing jet or an internally rotating jet flow. In the latter case, 
the expected observed period is typically P obs <10 days for 
massive quasars and even less for objects with significant 
bulk Lorentz factors, i.e., blazars. Given the range of values 
we have identified, this seems an unlikely model. 

In a single SMBH system, precession of the jet can arise 
if the accretion disk is tilted (or warped) with respect to the 
plane of symmetry (although the origin of the tilt would still 
need to be accounted for). Using data in the literature, Lu 


7 https://github.com/brandonckelly/carma.pack 


& Zhou (2005) find that the expected precession period for 
a jet with a single SMBH ( Lu &; Zhou|2005 ) is: 


log(r pre c/yr) ~ 0A8M abs + 0.14 log 


M 


10 8 M 0 


+ A 


where A = 15.18 — 18.86, depending on the values of con¬ 
stant parameters relating to the viscosity and angular mo¬ 
mentum of the SMBH system. For a representative quasar 
with M a bs — —25 with a 10 8 M© SMBH, this gives a pe¬ 
riod between 10 2 ' 2-6,9 years which is much longer than the 
observed periods reported here. We note as well that jets 
in non-blazar systems are predominantly associated with 
geometrically thick accretion flows in low-luminosity AGN 
(.Lboi ^ 10 42 erg s -1 ) whereas the candidate objects con¬ 
sidered here span the range 10 45 < Lboi < 10 48 erg s 1 . 
Reported jetted systems in the high luminosity range tend 
to be blazars and radio galaxies, e.g., Ghisellini et al. (2010), 


Sbarrato et al. (2014), but only a few of the candidates fit 


that category. It is still possible, though, that the central 
part of the accretion disk in the highest luminosity systems 
could be geometrically thick with a slim disk associated to 
accretion rates close to or above the Eddington limit, e.g. 
Sadowski & Narayan (2015) and references therein. 


Shorter periods are possible, however, for jet procession 
in a SMBH binary system. In this case, the jet precesses as 
a result of inner disk precession due to the tidal interaction 
of an inclined secondary SMBH. 


6.2.2 Warped accretion disk 


In addition to precessing any potential jet, a warped accre¬ 
tion disk might be responsible for the exhibited periodic¬ 
ity by obscuring the continuum-emitting region or at least 
modulating its luminosity as it precesses. Warped features 
are known to exist when the data is of sufficient quality to 


identify them, e.g. Herrnstein et al. (2005) for NGC 4258 


and may be common to most AGN, according to high reso¬ 


lution simulations of gas inflow in active galaxies (Hopkins 
Sz Quataert|2010 ). In general, the warping may be produced 
by a number of different mechanisms: the Bardeen-Peterson 
effect, a quadrupole potential or tidal interaction in a bi¬ 
nary system, self-gravity, angular momentum transport by 
viscous disk stresses, and radiation- and magnetic-driven in¬ 
stabilities. 

Simulations suggest that the Bardeen-Peterson effect 
may not typically occur in AGN disks ( Fragile et al.| 2007). 
Radiation- and magnetic-driven instabilities are also not fea¬ 
sible for AGNs as they predict behaviors inconsistent with 
observation ( Gaproni et al.|2006 ). Tremaine 8z Davis (2014) 
argue that, in the absence of any source of external torque, 
self-gravity is expected to play a prominent role in the dy¬ 
namics of AGN accretion disks. The precession rate of a self- 
gravitating warped disk around a SMBH is (Ulubay-Siddiki, 
Gerhard 8z Arnaboldi|2009 ): 


^ c( 


GM 2 d 

\ MbhGjj 


where r w is the radius of the warp, Mbh is the mass of the 
central object, M d is the mass of the accretion disk, and C is 
a constant of order unity accounting for the disk configura¬ 
tion (inclination, etc). For a fiducial value of Mbh — 10 8 M Q , 
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Tremaine & Davis (2014) find a warp radius of ~ 300 R g 
( R g = GMbh/c ) and, assuming Mdisk — 0.01 Mbh, the 
precession timescale is 51 years. This is only an order of 
magnitude larger than the typical periods reported here. The 
typical lifetime for a warp is talign — 1.3 x 10 5 (r w /300R g ) 4 
yr (Tremaine & Davis 2014), much shorter than the typi- 


6.2.5 Relativistic beaming 


cal AGN lifetime, and so it would also be a relatively rare 
phenomenon to encounter. 

A common explanation for a warped accretion disk is a 
binary system (MacFadyen & Milodavljevic 2008; Hayasaki 


et al.||2014 Tremaine &; Davis 2014), particularly in stellar 

systems where all components are resolvable. In the case of 
a SMBH binary, a warp will occur in the accretion disk of 
one of the SMBHs if its spin axis is misaligned with the 
orbital axis of the binary. It is also possible that there is a 
circumbinary accretion disk which is warped. 


6.2.3 Periodic accretion 


Farris et al. (2014) describe how periodic mass accretion 
rates in a SMBH binary system can give rise to an over dense 
lump in the inner circumbinary accretion disk (CAD) and 
that the observation of periodicity in quasar emission would 
be associated with this. Periodic accretion would not be ex¬ 
pected for a single SMBH system. [Hayasaki, Mineshige &; Ho| 
(2008) propose that small temperature changes associated 
with accretion variations would produce large-amplitude 
variations in the UV region, where the spectral energy distri¬ 
bution (SED) of the quasar decays exponentially, but would 
have little effect on the Rayleigh-Jeans part of the spectrum. 
However, [Rafikov (2013) argues that the SED of a circumbi¬ 
nary disk exhibits a power law segment with uF u oc v 12 / 7 
rather than v 4 ^ 3 . Accretion variations would therefore have 
a more noticeable effect at UV wavelengths in a binary sys¬ 
tem. We note as well that close pre-main-sequence binary 
stars show periodic changes in their observed optical lumi¬ 
nosity which is attributed to periodic accretion from the 


circumbinary disk, e.g., Jensen et al. (2007). 


D’Orazio et al. (2015) have also advanced an alternate sug¬ 


gestion for the periodicity of PG1302-102 - the emission seen 
is from a mini-disk around the secondary black hole which 
is in orbital motion around the system barycenter. Doppler 
boosting of the emission as the disk moves along the line-of- 
sight is sufficient to account for the periodic variation seen. 
With the right choice of parameters, the variability of PG 
1302-102 can explained by this if the broad lines do not orig¬ 
inate from any circumbinary disk. The model also predicts 
different amplitude variability at different frequencies, as¬ 
suming that the spectral slope is different, with no phase 
difference, i.e, it should track the optical variability. 


6.2.6 Quasi-periodic oscillation 

Quasi-periodic oscillations (QPOs) are a common phe¬ 
nomenon in the X-ray emission of stellar mass black hole 
binaries. Although the underlying physical cause is not un¬ 
derstood, high frequency QPOs are consistent with the dy¬ 
namical time of the system whilst low frequency ones are as¬ 
sociated with Lense-Thirring precession of a geometrically 
thick accretion flow near the primary black hole, e.g., |In-| 
gram , Done V Fragile| |2009). If the frequencies involved are 
scaled up by mass (many black hole timescales depend ap¬ 
proximately on the inverse of the black hole mass), they are 
roughly consistent with the periodic behavior reported in 
this work. For example, the microquasar GRS 1915+105 has 
a mass of ~ 12M© and exhibits QPOs at ~1 Hz (e 


Yan 


et al. 2013 and references therein) with slowly varying fre¬ 


quency. Scaling this to the observed mass range of PG 1302- 
102 (10 8 ' 3_9 ' 4 Mq) gives QPOs with a period between ^200 
and 2400 days (the restframe period is ~1440 days). King et 
al. (2013) reported the first AGN analog of a low-frequency 


QPO in the 15 GHz light curve of the blazar J1359+4011 
with a timescale varying between 120 and 150 days over a ~4 
year timespan. Stellar black hole binary QPOs are predicted 
to show a decreasing rms amplitude with longer wavelength 
(Veledina V Po utanen|2015 ). The degree and angle of linear 
polarization of the precessing accretion flow are also pre¬ 


dicted to modulate on the QPO frequency (Ingram et al 


2015). 


6.2.4 Accretion disk gap 

From hydrodynamical simulations of SBMH binaries, 


D’Orazio et al. (2015) have proposed that the detected pe¬ 
riodicity in PG1302-102 could be associated with a lopsided 
central cavity in the CAD for q > 0.3, where q is the ra¬ 
tio of the two SMBH masses: q = M 2 /M 1 . Significantly, in 
this scenario, the strongest period detected is a factor of 
X = 3 — 8 times greater than the true period of the bi¬ 
nary and so the true binary separation is a factor of x _2/ ^ 
smaller. This model predicts additional smaller amplitude 
periodic variability on timescales of tu n and 0.5 tbin ( hin is 
the true binary periodicity), and periodic spectral variations 
in broad line widths and narrow line offsets. There are also 
measurable relativistic effects on the Fe K a line and an ob¬ 
servable decrement in the spectral energy distribution of the 
source at a wavelength determined by the width of the gap, 
e.g., Giiltekin & Miller ( |2012| . 


6.3 Theoretical detection rates 

Given our large data set of 250,000 quasars, it is worthwhile 
considering how many SMBH binary systems we could ex¬ 


pect to find. Haiman, Kocsis & Menou (2009) used simple 


disk models for circumbinary gas and the binary-disk inter¬ 
action to determine the number of SMBH binaries that may 
be expected in a variety of surveys, assuming that such ob¬ 
jects are in the final gravitational wave-dominated phase of 
coalescence (which equates t o separations less t han ^0.01 


pc for a 10 8 M o system). Volonteri, Miller & Dotti (2009) 
combined this with merger tree assembly models to similarly 
predict the number of expected SMBH binaries at wider sep¬ 
arations where spectral line variations may be seen (equating 
to separations greater than ~ 0.2 pc for a 1O 8 M0 SMBH bi¬ 
nary). The latter shows that in a sample of 10000 quasars 
with z < 0.7, there should be ~ 10 such objects and this 


© 2011 RAS, MNRAS 000, [I]-?? 

































































26 M. J. Graham et al. 


number increases by a factor of ~5—10 for z < 1. We note 
that our sample has ~ 75000 quasars with z < 1. 

Using these approaches, we can estimate our predicted 
sample size. Assuming a limiting magnitude of V ~ 20, 
a detectable range of orbital periods from 20-300 weeks 
(thus spanning both gravitational-wave and gas-dominated 
regimes), a survey sky coverage of 2 i r steradians (there are 
far fewer known spectroscopically confirmed quasars in the 
regions where CRTS does not overlap with SDSS), and a 
redshift range of 0.0 - 4.5, we should identify ~ 450 SMBH 
binaries. Our detection rate of 10 -4 is therefore conservative 
compared to the predicted rate of 2.3 x 10 -3 . 

Virial black hole mass estimates (Shen et al. 2011) are 
available for 88,000 quasars in our sample (of which 23 per 
cent have z > 2). If we assume that each of these is a SMBH 
binary with a separation of 0.01 pc then the CRTS temporal 
baseline is sufficient to detect 1.5 cycles or more of periodic 
behavior in 63 per cent of them (including 55 per cent of 
the z > 2 population). Our search should therefore be sen¬ 
sitive to a large fraction of the close SMBH population. We 
note, however, these theoretical arguments are still subject 
to considerable uncertainties, for example, if the final par¬ 
sec problem cannot be resolved then there will not be any 
binaries in the ~0.01 pc regime. 


Haiman, Kocsis & Menou (2009) also predict that in the 


gravitational wave (GW)-dominated regime of the merger 
process, the number of expected binaries is: 


/10 7 yr\ 

torb 

V tQ ) 

50.2 week 


8/3 


Mr 


- 5/3 


Qs 


where t or b is the restframe period, tQ is the timescale, M 7 
is the black hole mass in units of 10 7 M© and q s is the 
black hole mass ratio. Fig. [9] shows the observed distribu¬ 
tions for the binary candidates with log (Mbh /Mq) < 9 and 
\og(MsH/M q ) > 9 together with n va r oc t 3 (! 3 relationships 
fit to these (log (Mbh/Mq) ~ 9 is the median value for the 
data set and provides a natural division point). The results 
suggest that there is a transition timescale below which ob¬ 
jects follow the power law expected from orbital decay driven 
by GWs. This timescale also anticorrelates with object mass 
such that higher mass objects (here log (Mbh /Mq) > 9) 
change behavior at about t or b ~ 750 days and lower mass 
objects at t or b ~ 1000 days. 

At longer timescales (and lower masses), the number 
of objects is expected to follow a power law distribution 
n var oc £ a , where the scaling index a is dependent on the 
physics of the circumbinary accretion disk and viscous or¬ 
bital decay. However, since one of our selection criteria is 
that the temporal coverage of a light curve should cover 
at least 1.5 cycles of the periodic signal in the observed 
frame, we are incomplete at restframe timescales larger than 
the transition timescale. We can therefore not say anything 
about the true value of a and the physics of the gas-driven 
phase from our sample due to insufficient coverage. Never¬ 
theless the statistical detection of GW-driven binaries can 
also be seen to confirm that circumbinary gas is present at 
small orbital radii and is being perturbed by the black holes 
(Haiman, K ocsis &; Menou|[2 009). 


Figure 9. The frequency of supermassive black hole binaries as 
a function of restframe period. This shows the distribution for 
objects with log (Mbh /Mq) < 9 (blue) and log (Mbh /Mq) > 
9 (green). Theoretical curves are also shown for a n va r oc £ 3 ^ 3 
relationship for the two mass distributions: \og(MBH /Mq) < 9 
(black) and log (Mbh/Mq) > 9 (red). 



7 GRAVITATIONAL WAVE DETECTION 


Given that some of the candidates appear to be in the 
gravitational-wave driven phase of merging, it is worth de¬ 
termining whether any of them might be resolvable as indi¬ 
vidual sources for current or upcoming nanohertz GW de¬ 
tectors (pulsar timing arrays (PTA)) rather than just be 
contributors to the stochastic GW background. For each 
candidate, we have determined the maximum signal that 
could be detected (see Table [2| . The GW frequency of a bi¬ 
nary with a circular orbit is few = 2 /t per , where t per is the 
orbital period. The intrinsic GW strain amplitude is: 


where M c = (MiM 2 ) 3 / 5 (Mi+M 2 ) _1//5 is the chirp mass and 
dh is the comoving distance to the source. The maximum 
induced timing residual amplitude is then h s /(‘lit few)- 

We have also calculated the inspiral time for each can¬ 
didate (see Table [ 2 ]) to check that there is no significant fre¬ 
quency evolution over the lifetime of a detection experiment, 
i.e., t eX pt tinspi where ti nsp = 5(1 + q) 2 a A /2qR% and 
R s — 2 GM/c 2 . The median inspiral time is ~8,000 years, 
although there are four candidates predicted to merge within 
the next century and , with decadal baselines, we should be 
able to detect period changes in these objects. 

Using PTA noise estimates for Nanograv and the Parkes 
PTA (Arzoumanian et al. 2014 Zhu et al. 2014), we calcu¬ 
late that none of the candidates will be resolvable as sources 
by the current generation of experiments (see Fig. 10). How¬ 
ever, some of the sources will be viable by the next genera¬ 
tion of detectors, e.g., SKA. 


8 CONCLUSIONS 

We have detected strong periodicity in the optical photom¬ 
etry of 111 quasars which we ascribe to the presence of a 
close SMBH binary in these systems. The Keplerian nature 
of the signal suggests that it is kinematic in origin and may 
be produced by a precessing jet, warped accretion disk or 
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Figure 10. The expected GW strain amplitude for the candidate 
sources. The error bars indicate the range of strain amplitudes for 
q = 0.05 to q = 1. The blue line indicates the all-sky averaged 
noise limit for the Nanograv experiment and the red line that of 
the most sensitive detection limit for the Parkes PTA. 



log(Frequency (Hz)) 


periodic accretion driven by the SMBH binary. We note our 
detection methodology may not be particularly sensitive to 
other types of periodic behavior exhibited by SMBH binary 
candidates, such as OJ 287. This would then suggest that 
our sample of objects represents a small fraction of a much 
larger close binary SMBH population which is yet to be 
detected. We are therefore planning further studies of the 
(quasi-)periodic quasar population with less stringent con¬ 
straints than those applied in this analysis. This will also 
address our incomplete coverage of the gas-driven merger 
regime. 

Further observations of our sample can test the differ¬ 
ent interpretations, as well as the primary SMBH binary 
attribution. Shen V Loeb (2010) suggest that reverbera¬ 
tion mapping is a particularly useful diagnostic since the 
behaviour of emission line response to continuum variations 
is expected to be different for alternate explanations. We 
plan to see whether stacked reverberation employing time 
series and pairs of multiepoch spectra (Fine et al. 2013) 
is sufficient. Obviously continued monitoring by CRTS and 
other synoptic surveys will track future cycles of periodic¬ 
ity and historical photometric data from photographic plate 
collections, such as DASCfJ] may provide more data for 
previous cycles. For example, DASCH data for PG 1302-102 
is consistent with the reported period (J. Grindlay, private 
communication). With such decadal baselines, any change 
in the period of the system expected from general relativity 
may be detectable. 

Future spectroscopic observations can also test whether 
there is any spectral variation in the sample consistent with 
binary orbital timescales, although this is not necessarily 
present in such close systems. The theoretical expectation 
for pairs at the close separations that we are probing is that 
they should not show any such effect - the size of any broad 
line region (BLR) dwarfs the orbital dimensions of the bi- 


8 http://dasch.rc.fas.harvard.edu/project.php 
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nary and at separations closer than 0.03 pc in some sce¬ 
narios, the BLR is truncated or destroyed ( Ju et ah||2013 ). 
However, D’Orazio et al. (2015) suggest that there may be 
broadened emission due to recombination in orbiting cir- 
cumbinary gas which would show as periodic variation in 
the FWHM of associated lines and smaller shifts in their 
centroids. Continued spectral monitoring, particularly co¬ 
incident with an extremum in the photometric light curve, 
will help to confirm the origin of any spectral variability de¬ 
tected. With a binary, the variation should follow a regular 
evolution whereas a bright spot (non-axisymmetric pertur¬ 
bation in the BLR emissivity) in the BLR, say, should just be 
a transient feature. One caveat with any detection of spec¬ 
tral variability is that asymmetric reverberation can act as a 
major source of confusion noise in multiepoch spectroscopic 
data (Barth et al.|2015 ). 

Multiwavelength observations, particularly of the ra¬ 
dio quiet objects, will also provide more information about 
them. In particular, the specific predictions that |D’Orazio et] 


al. (2015) make regarding the Fe Ka line, broad line widths 


and narrow line offsets for the accretion disc cavity model 
can be tested. Finally, we note that these objects are strong 
candidates for any gravitational wave experiment sensitive 
to nanohertz frequency waves such as those using pulsar tim¬ 
ing arrays as well as any future space-borne gravitational 
wave detection mission. 
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